// Fracture width
d = 0.01;

// Domain dimensions
Lx = 1;
Ly = 0.5;
Sf = 0.1;

// Mesh size
h  = 0.05;
hf = d/2.0;

// Define the reference points
A = newp; Point(A) = {-Lx,-Ly,0.0,h};
B = newp; Point(B) = { Lx,-Ly,0.0,h};
C = newp; Point(C) = { Lx, Ly,0.0,h};
D = newp; Point(D) = {-Lx, Ly,0.0,h};

E = newp; Point(E) = {-d/2.0,-Ly,0.0,hf};
F = newp; Point(F) = { d/2.0,-Ly,0.0,hf};
G = newp; Point(G) = { d/2.0, Ly,0.0,hf};
H = newp; Point(H) = {-d/2.0, Ly,0.0,hf};

I = newp; Point(I) = {-d/2.0,-(Ly-Sf),0.0,hf};
L = newp; Point(L) = { d/2.0,-(Ly-Sf),0.0,hf};
M = newp; Point(M) = { d/2.0, (Ly-Sf),0.0,hf};
N = newp; Point(N) = {-d/2.0, (Ly-Sf),0.0,hf};

// Define the five subdomains
b1 = newc; Line(b1) = {A,E}; 
b2 = newc; Line(b2) = {F,B}; 
t1 = newc; Line(t1) = {H,D}; 
t2 = newc; Line(t2) = {C,G}; 
rr = newc; Line(rr) = {B,C}; 
ll = newc; Line(ll) = {D,A}; 

bf = newc; Line(bf) = {E,F}; 
tf = newc; Line(tf) = {G,H}; 

sb1 = newc; Line(sb1) = {E,I}; 
st1 = newc; Line(st1) = {N,H}; 
fl  = newc; Line(fl)  = {I,N}; 

sb2 = newc; Line(sb2) = {L,F}; 
st2 = newc; Line(st2) = {G,M}; 
fr  = newc; Line(fr)  = {M,L}; 

fb  = newc; Line(fb)  = {I,L}; 
ft  = newc; Line(ft)  = {M,N}; 

bO1 = newreg; Line Loop(bO1) = {b1,sb1,fl,st1,t1,ll}; 
bO2 = newreg; Line Loop(bO2) = {b2,rr,t2,st2,fr,sb2}; 
bSb = newreg; Line Loop(bSb) = {bf,-sb2,-fb,-sb1}; 
bSt = newreg; Line Loop(bSt) = {-ft,-st2,tf,-st1}; 
bFf = newreg; Line Loop(bFf) = {fb,-fr,ft,-fl}; 

Omega1 = news; Plane Surface(Omega1) = {bO1};
Omega2 = news; Plane Surface(Omega2) = {bO2};
SideB  = news; Plane Surface(SideB)  = {bSb};
SideT  = news; Plane Surface(SideT)  = {bSt};
Fract  = news; Plane Surface(Fract)  = {bFf};

// Display boundary labels
//View "Boundary labels" {
//  T2(10, 20, 0){ StrCat( 
//    "Hole labels:    ",Sprintf("%.0f, %.0f, %.0f, %.0f",a1,a2,a3,a4)) };
//  T2(10, 40, 0){ StrCat( 
//    "Inflow label:   ",Sprintf("%.0f",lef)) };
//  T2(10, 60, 0){ StrCat( 
//    "Outflow labels: ",Sprintf("%.0f, %.0f, %.0f",bot,rig,top)) };
//};

